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1.  Introduction 

This  paper  examines  in  some  detail  the  distribution  of  the  total  number  of 
cases  in  an  epidemic  of  the  general  stochastic  type  for  a  closed  population.  The 
assumed  model  is  that  of  Bartlett  [2]  and  McKendrick  [11]  which  Bailey  [1] 
used  to  study  the  stochastic  analogue  of  the  deterministic  threshold  theorem 
(Kermack  and  McKendrick  [10],  D.  G.  Kendall  [9]).  Bailey  obtained  recurrence 
relations  from  which  the  required  probabilities  were  computed  numerically. 
His  calculations  revealed  a  gradual  transition  from  ./-shaped  distributions  con¬ 
taining  only  small  epidemics  for  population  sizes  below  the  threshold,  to  U- 
shaped  distributions  containing  either  large  or  small  epidemics  but  practically 
no  epidemics  of  intermediate  size  when  the  threshold  is  exceeded.  There  is  also 
an  interesting  transitional  form  of  distribution  near  the  threshold  value. 

In  an  attempt  to  understand  what  motivates  an  epidemic  to  behave  in  this 
way,  Whittle  [13]  and  Kendall  [9]  constructed  different  models  approximating 
to  the  one  used  by  Bailey  but  easier  to  handle  analytically.  Both  explained 
Bailey’s  results  in  terms  of  an  initial  birth  and  death  process  where  extinction 
is  certain  in  the  first  case  and  not  certain  in  the  second.  This  work  is  summarized, 
with  additional  references,  in  the  book  by  Bailey  [2].  In  a  paper  presented  at 
this  Symposium,  Gani  [7]  develops  some  recent  work  by  Siskind  [12]  and 
himself  [6]  on  a  method  of  obtaining  time  dependent  solutions  of  the  epidemic 
equations.  For  the  limiting  case  considered  here  he  shows  how  the  probabilities 
can  be  computed  by  successive  multiplication  of  matrices. 

My  main  object  is  to  arrive  at  approximate  formulae  for  the  distribution  of 
total  epidemic  size  which  are  appropriate  for  large  populations.  The  approach 
differs  from  that  of  most  other  investigations  in  that  the  backward  equations 
of  the  process  are  used.  (See,  however,  Bartlett  [3].)  I  also  find  it  convenient  to 
work  in  terms  of  the  number  remaining  uninfected,  rather  than  the  total  number 
of  new  cases.  The  technique  by  which  the  approximations  are  obtained  was 
previously  used  by  me  in  an  entirely  different  context  (Daniels  [4],  [5]).  As 
presented  here  it  should  not  be  regarded  as  rigorously  establishing  the  approx¬ 
imations,  though  I  have  no  doubt  that  the  results  are  correct  and  numerical 
comparisons  bear  this  out. 
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2.  The  deterministic  model 


Suppose  that  at  time  t  there  are  x  susceptibles,  y  infectives  and  z  recovered 
or  dead  in  the  population.  Initially  we  shall  take  x  =  £,  y  =  77,  z  =  0,  so  that 
x-\-y  +  z  =  ^-\-ij.  The  deterministic  epidemic  has  been  fully  studied  by 
Kendall  [9],  but  we  examine  it  briefly  for  the  sake  of  some  results  referred  to 
later.  The  deterministic  equations  are,  with  a  suitable  time  scale, 


(2.1) 


dx  dy  dz 

di  “  ~xy’  Tt=xy~n’  di  = 


where  p  is  called  the  relative  removal  rate  by  Bailey,  and  the  threshold  by  Kendall. 
Then  dx/dz  =  —x/p  and  x  =  £  exp  (—z/p)  for  all  t.  At  the  end  of  the  epidemic 
y  =  0,  2  =  £  +  77  —  x,  and  the  number  x  of  individuals  remaining  uninfected 
satisfies  the  equation 

(2.2)  x  exp  (~x/p)  =  £  exp  [-  (£  +  7 ;)/p]. 

We  suppose  that  £  and  p  are  large  and  77/p  is  small.  There  are  two  values  of  x 
satisfying  (2.2)  near  the  respective  roots  £,  £'  of  x  exp  (—x/p)  =  £  exp  (— £/p). 
The  only  relevant  root  is  the  one  less  than  £.  When  £  <  p,  this  root  is  near  £ 
and  (2.2)  gives  as  a  first  approximation 

(2.3)  x  =  £  —  tj£/(p  —  £). 

On  the  other  hand,  when  £  >  p  the  required  root  is  near  £'  <  p  and  approx¬ 
imately, 

(2-4)  x  —  g  —  77 £'/(P  -  £')• 

These  hold  provided  £  is  not  near  p,  though  it  is  easy  to  get  an  approximate 
transitional  form  for  £  ~  p.  Together,  they  constitute  the  deterministic  thresh¬ 
old  theorem  in  terms  of  the  number  remaining  uninfected.  Notice  that  when 
£  ^ >  p  then  £'/p  is  small  and 

(2.5)  x  ~  £'  ~  £e-t//>. 


3.  The  stochastic  model 

In  the  continuous  time  model  considered  by  Bartlett  and  McKendrick  the 
transitions  in  (t,  t  +  8t)  are  (x,  y)  — >  (x  —  1,  y  +  1)  with  probability  xyU  + 
o(8t)  and  (x,  y)  — >  (x,  y  —  1)  with  probability  py8t  +  o(8t).  As  we  are  concerned 
only  with  the  final  distribution  of  x ,  it  is  simpler  to  work  with  the  random 
walk  of  transitions  in  the  x,  y  plane  such  that 

(3  P{(*,  y)  -» (x  ~  1,  V  +  1)}  =  */(p  +  x), 

P{(x,  y)  -»  (x,  y  -  1)}  =  p/(p  +  x), 
where  initially  x  =  £,  y  =  77  and  absorption  occurs  on  y  =  0. 
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An  alternative  formulation  of  the  random  walk  in  terms  of  the  numbers  of 
new  cases  w  =  £  —  x  and  removals  z  is  of  some  interest.  This  is 

(3  2)  P^W'  ^  ^  +  =  ^  ”  ™)/(p  +  £  “  w)> 

P{(M,2)->(ld,2+l)}=p/(p  +  ^-w). 

The  problem  can  then  be  described  in  terms  of  a  game  involving  a  mixture  of 
sampling  with  and  without  replacement.  A  box  contains  £  real  pennies  and  p 
false  ones.  The  player  starts  with  a  capital  of  rj  pennies  and  the  price  of  a  draw 
is  one  penny.  If  he  draws  a  false  penny  he  replaces  it  in  the  box.  If  he  draws  a 
real  penny  he  keeps  it  and  is  allowed  a  further  trial  free.  The  game  stops  when 
the  player  is  ruined  (w  +  -q  =  z)  or  when  he  has  drawn  all  the  real  pennies 
(w  =  £). 

Let  p(x |£,  77)  be  the  probability  that  there  are  ultimately  x  uninfected  individ¬ 
uals,  when  initially  there  were  £  susceptibles  and  7 7  infectives.  The  backward 
equations  for  p  are 

(3.3)  £p(z|£  —  1,  m  +  1)  +  pp(x |£,  77  -  1)  —  (p  +  £)p(z|£,  17)  =  0, 

£  >  x,  77  ^  1,  and 

(3.4)  pp(x \x,  77-I )  -  (p  +  x)p(x \x,  7?)  =  0 
with  the  condition 

(3.5)  p{x |£,  0)  =  5(£  -  x), 
where  5(£  —  x)  =  0,  £  5^  x  and  5(0)  =  1.  From  (3.4), 

(3.6)  p(x \x,  77)  =  [p/(p  +  x)]p(x \x,  77  -  1)  =  •  •  •  =  [p/(p  +  *)]’. 

Our  method  of  attack  depends  on  the  fact  that 

<«  (.iJC-rfr-.r’ 

is  a  solution  of  (3.3)  for  arbitrary  s,  and  satisfies  (3.6)  when  s  =  0.  We  try  to 
build  up  a  solution  of  the  form 


(3.8)  P(x\k,  ,)  =  %  A.  (x  l  )  (p-Tf+1)S"+’ 

satisfying  the  required  conditions.  If  A0  =  1,  condition  (3.6)  is  satisfied,  and 
from  (3.5)  we  must  have 

(3.9)  »({ -X)  =  l :  A,  (  *  )  (  p  Y-\ 

«  =  0  \Z  +  S/  \p  +  x  +  s/ 

The  coefficients  Aa  can  be  determined  recursively  and  hence  p (z|£,  77)  can  be 
found,  provided  Aa  is  independent  of  £.  That  this  is  so  becomes  clear  if  we  write 


(3.10) 
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and  use  the  fact  that 


(3.11) 
Then, 

(3.12) 

(3.13) 


The  coefficients  H,  depend  only  on  x  and  p. 

As  a  method  of  computing  the  probabilities,  these  equations  are  if  anything 
less  convenient  than  the  original  equations  (3.3),  (3.4),  or  the  corresponding 
forward  equations,  when  £  is  large  and  £  >  p.  Their  value  lies  in  the  fact  that  a 
technique  is  available  for  obtaining  an  asymptotic  approximation  when  £  is 
large. 


4.  Some  exact  results 


The  problem  is  substantially  simplified  by  using  the  following  result.  Since 
the  left  side  of  (3.13)  is  zero  except  when  £  =  x,  it  can  equally  well  be  put  in 
the  form 


(4-1) 


*(£  -  *)  =  L 
8  =  0 


p  +  x  y-* 
p  4*  X  +  s) 


If  we  then  write  H„  =  H,(x,  p)  and  p(x\£,  17)  =  p(x |£,  rj,  p)  to  show  their  de¬ 
pendence  on  p,  it  appears  that 


(4.2)  H„(x,  p)  =  H,{ 0,  p  +  x ), 

and  from  (3.13)  we  get  the  exact  relation 


(4.3)  p(x  |£,  r j,  p)  =  ,(^)  P(°li  ~  X>V  >P  +  x)- 

(H.  F.  Downton  has  succeeded  in  deriving  (4.3)  by  a  direct  combinatorial 
argument.)  This  enables  everything  we  want  to  be  deduced  from  a  knowledge 
only  of  the  behavior  of  p(0|£,  7?,  p)  which  we  now  study.  The  equations  are, 
with  Hg  for  Hs( 0,  p), 


(4.4) 

p(0|$,  v ,  p)  =  E  (- 

-)'H( 

8=0 

(4.5) 

m  =  l  (- 
8  *0 

- )'H , 

Taking  £  =  0,  1,  2,  •  •  •  in  (4.5)  we  get  a  set  of  equations  for  Ha  which  can  be 
solved  in  determinantal  form  as 
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si 


(4.6)  h.  = 


Co 

4 

2! 

4 

3! 


cS"1 


Cl 

1! 

Ci 

2! 


csr2 


o  oo 

1  o  o 

c2  10 


Cl 


(«-!)!  (s  —  2)!  (8-3)! 


Co 
s ! 


c-r1 


cl  2 


Cs— 2 


Cs  —  2 


Cg— 1 


(8-1)!  (8-2)!  2! 

where  c»  =  p/(p  +  s)  and  =  1.  Substitution  in  (4.4)  gives  an  explicit  solu¬ 
tion  for  p(0|£,  t/,  p). 

6.  Approximations  below  the  threshold 

The  form  of  the  solution  (4.6)  bears  an  unexpected  and,  I  suspect,  fortuitous 
resemblance  to  expressions  for  the  probabilities  associated  with  the  Poisson 
process  with  a  curved  absorbing  boundary  and  the  related  Kolmogorov-Smirnov 
test  [5].  For  example,  it  follows  from  the  work  on  that  problem  that  H„  can  also 
be  expressed  in  the  form  Hs  =  cs~s  C,( 0),  where 


(5.1) 


C8(z)  —  s!  J™  dzi  f*dz2  •  •  •  f°'* 


dz , 


is  a  so  called  Gontcharoff  polynomial  [8]  whose  jth  derivative  vanishes  at 
z  =  Cj.  A  good  deal  can  be  discovered  about  the  asymptotic  behavior  of  these 
expressions  by  using  a  technique  originally  developed  for  a  related  problem  [4]. 
But  there  is  an  essential  difference  here  which  complicates  matters.  In  the 
applications  mentioned  c*  is  an  increasing  positive  sequence  and  this  ensures 
that  C«(0)  is  always  positive.  In  the  present  problem  c«  is  a  decreasing  positive 
sequence  and  beyond  a  certain  value  of  s,  Hs  begins  to  oscillate  with  increasing 
amplitude  and  alternating  sign.  Nevertheless,  we  shall  use  a  modified  version 
of  the  same  technique  to  study  the  asymptotic  form  of  Hs  and,  hence,  deduce 
that  of  p(0|£,  77,  p)  when  £  and  p  are  large.  The  presentation  is  heuristic  and  to 
some  extent  incomplete.  A  rigorous  development  must  depend  on  a  more  ex¬ 
tensive  study  of  the  asymptotic  properties  of  Gontcharoff  polynomials  which 
it  is  hoped  to  publish  later. 

In  (4.5)  replace  £  by  m,  multiply  it  by  (  —  X)m  and  sum  from  m  =  0  to  £. 

After  a  little  manipulation  the  result  is 


(5.2) 


1  =  £  (Xc.)*(l  ~  Ac.)«  *,  c ,  =  p/(p  +  s), 

8=0  \S/ 
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where  X  is  an  arbitrary  parameter.  The  technique  is  to  look  at  the  behavior  of 
(5.3)  r.(X)  =  (*)  (Xc.)-(l  -  Xc.)f- 


for  large  £.  It  can  be  shown  to  have  limiting  forms  analogous  to  the  normal 
and  Poisson  limits  of  the  binomial  distribution,  with  a  peak  at  the  unique  root 
s0  of  s  =  X£c„  By  varying  X  we  can  scan  Hs  with  this  “window”  and  deduce  its 
asymptotic  behavior.  The  range  0  <  X  <  1  ensures  that  Ts(\)  is  positive,  but 
X  can  be  allowed  to  exceed  unity  provided  Ts{\)  remains  positive  near  the 
root  So. 

We  shall  consider  only  the  normal  limiting  form.  Assume  £  and  p  are  large 
and  write  z  —  s/£,  dz  =  l/£,  c(z)  =  ct,  Ts(\)  =  T(z,  \)dz.  If  neither  z  nor  1  —  z 
is  small,  Stirling’s  approximation  leads  to 


(5.4)  T(z,  X)  ~  [; 


£  i1'2  rxc(z)"p2  ri  -  xcCzn^1-*) 

J  L  i  - 


2x2(1  —  z)  J  L  2  J  L  1  -  2  J 

[z  -  Xc(z)]  1  f  f 


j  ,  [e  -  Xc(z)]y»-2> 
1-2  J 


r  j  yn  \  t[l  _  Xc(2)]!  ,  involving^ 

=  L^i) J  + 0h;f ”  } 


The  maximum  of  T(z,  X)  is  at  the  unique  root  z0  of  z  —  \c(z)  =  0,  and  it  is 
shown  in  [4]  that  there  are  no  other  maxima.  Near  z0  we  can  write  z  —  Xc(z )  4= 
(z  —  z0)[l  —  Xc'(zo)],  and  because  z  —  z0  is  0(£1/2)  over  its  effective  range,  we 
get  the  normal  approximation 


(5.5) 


rfex)~[2^a-Je*p{- 


£[1  -  Xc’(z0)]2 
2z0(l  —  z0) 


where 


(5.6)  Zo  =  Xc(z0). 

If  it  can  now  be  assumed  that  H(z)  =  H(Z  varies  slowly  with  z  near  z0,  then 
(5.2)  approximates  to 

(5.7)  1  ~  P  H(z)  T(z,  X)dz  ~  H(z„)  JQ'  T(z,  \)dz 

~ff(so)/[l  -  Xc'(zo)]  =  H(za)/[1  -  zoc'(zo)/c(z0)]. 

Hence, 

(5.8)  H(z)  ~  1  —  zc'(z)/c(z), 

provided  z  is  a  possible  root  of  (5.6)  at  which  (5.5)  holds.  In  terms  of  s,  we  have 
c'(z)/c(z)  =  —  £/(p  +  s),  and  we  arrive  at  the  approximation 

(5.9)  Ht  ~  1  +  s/(p  +  s)  —  2  —  p/(p  +  s). 

The  result  is,  of  course,  suggested  rather  than  established  by  this  kind  of 
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reasoning.  It  depends  on  the  assumption  that  H{z)  varies  slowly,  and  we  would 
expect  it  to  fail  for  large  values  of  s  corresponding  to  roots  for  which  X  causes 
Ts(k)  to  oscillate.  Also,  small  values  of  s  have  been  excluded  by  the  argument. 
Nevertheless,  calculations  show  that  the  approximation  is  good  for  values  of 
s  up  to  about  p  (see  table  I).  It  should  therefore  provide  an  approximation  for 
p(x  |£,  rj,  p)  when  £  is  less  than  the  threshold  p. 


TABLE  I 

Exact  and  Approximate  Values  H,  for  p  =  20 
Variation  of  K,  with  s 


s 

H. 

1  +  s/(p  +  s) 

K. 

0 

1.0000 

1.0000 

1.0000 

1 

1.0500 

1.0476 

-21.0000 

2 

1.0948 

1.0909 

109.4762 

3 

1.1352 

1.1304 

-336.3698 

4 

1.1721 

1.1667 

732.5702 

5 

1.2059 

1.2000 

-1234.8024 

6 

1.2369 

1.2308 

1696.7042 

7 

1.2655 

1.2593 

-1966.9936 

8 

1.2921 

1.2857 

1971.5763 

9 

1.3168 

1.3103 

-1740.1758 

10 

1.3397 

1.3333 

1371.8941 

11 

1.3612 

1.3548 

-977.0948 

12 

1.3813 

1.3750 

634.5692 

13 

1.4002 

1.3939 

-378.7168 

14 

1.4179 

1.4118 

209.0655 

15 

1.4347 

1.4286 

-107.3579 

16 

1.4504 

1.4444 

51.5285 

17 

1.4655 

1.4595 

-23.2195 

18 

1.4792 

1.4737 

9.8550 

19 

1.4938 

1.4872 

-3.9585 

20 

1.5028 

1.5000 

1.5028 

21 

1.5268 

1.5122 

-0.5480 

22 

1.4979 

1.5238 

0.1840 

23 

1.6522 

1.5349 

-0.0664 

24 

1.1344 

1.5455 

0.0143 

25 

3.1684 

1.5556 

-0.0120 

26 

-4.8683 

1.5652 

-0.0053 

27 

28.2780 

1.5745 

-0.0086 

28 

-112.890 

1.5833 

-0.0091 

29 

508.531 

1.5918 

-0.0106 

30 

-2314.03 

1.6000 

-0.0121 

31 

10901.2 

1.6079 

-0.0137 

32 

-52810.9 

1.6154 

-0.0155 

33 

263176 

1.6226 

-0.0175 

34 

-1347561 

1.6296 

-0.0197 

35 

7084284 

1.6364 

-0.0221 

Let  us  substitute  (5.8)  in  the  right  side  of  (4.5)  to  see  how  nearly  it  is  satis¬ 
fied  for  £  >  0.  (It  is  exact  for  £  =  0).  We  have 
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»■»>  a  <->•»•  o  te)'  ~k  <->■  (!)  te)'  -  ten 


(t-  l)l/o 


(2w*_1  —  w*/£) 


(1  —  e~ulp)t  e~udu 

-<rh 

—  (l/^  +  £/p)w2*/p*  +  •  •  •]  dw. 

The  expansion  within  the  integrand  is  convergent  for  small  enough  w,  and 
termwise  integration  yields  an  asymptotic  expansion  in  powers  of  p-1  for  fixed 
£.  The  term  in  p~f  vanishes  and  the  leading  term  is 

(5.11)  -*(2*)!/p*«(*-  1)!. 

For  large  p  this  is  small  even  at  £  =  1.  As  £  increases  it  becomes  approximately 

(5.12)  —  [pe2/8(2)1/2](4£/pe){+2 


which  decreases  to  a  minimum  at  about  £  =  p/4.  It  does  not  become  appreciable 
again  until  £  approaches  pe/4  after  which  it  rapidly  becomes  large.  So  at  least 
for  £  <  pe/4  one  can  with  some  confidence  insert  (5.9/  into  (4.4)  and  obtain  in 
the  same  way, 

a 

(5.13)  p(0fc  o,  p)  ~  ^-+  *  _  ^  Jo  [2 _„«+,/({  + ,)] 


(1  —  e~u,p)*e~u  du 


t?(2£  +  ^  ~  1)! 
(£  +  v)  1  Pf 


+  0(p-«-i). 


Then  from  (4.3)  we  get  for  values  of  £  below  the  threshold  p, 


(5.14) 


p(*l£»  v,  p) 


VP 


Z— -c+1? 


(p  +  z)2*-2^” 


(!) 


(2£  —  2a:  +  j?  —  1)! 


VP 


(p  4-  £  -  w)2^’ 


© 


(£  —  £  +  77)! 

(2^  +  7?  -  1)! 

(w  +  17) ! 


in  terms  of  the  number  w  of  new  cases.  Since  £  is  large,  a  further  approximation 
leads  to  the  result 


(5.15) 


qp«+*£”  (2w  +  77  —  1)! 

^  (p  +  £)2“’+’!  w>!(ip  +  77)! 


which  is  the  solution  for  the  birth  and  death  process  proposed  by  Bartlett  and 
exploited  by  Whittle  and  Kendall,  having  birth  rate  £,  death  rate  p,  and  whose 
mean  is  given  by  the  deterministic  approximation  (2.3). 
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6.  Approximations  above  the  threshold 


To  investigate  the  range  of  £  above  the  threshold  p,  it  is  necessary  to  study 
the  behavior  of  Hs  at  values  of  s  beyond  those  which  can  be  reached  by  the 
previous  method.  It  has  been  mentioned  that  Ha  begins  to  oscillate  violently 
when  s  exceeds  a  certain  value.  This  suggests  that  the  following  transformation 
of  (5.2)  will  be  useful  for  such  values  of  s.  Write  X  =  —v/(l  —  v ),  and 

(6.1)  H.  =  (-)•[(  1  -  csy/ct]Ks  =  ( -)•(?/ PyK. . 


Then,  (5.2)  becomes 

(6.2)  (1  -»)«  =  £  K.  (f)  [,(1  -  0]*[1  -  *(1  -  *.)]«-. 

We  could  now  try  using  the  previous  technique  with  Tt(\)  replaced  by 

(6.3)  V. w  =  (9  0(1  -  c.)]>[l  -  f(1  -  c.)]«-, 
which  has  peaks  at  the  roots  of 

(6.4)  s  =  »>£(1  —  c8)  =  vs£/(p  +  s). 

The  lower  root  s  =  0  is  irrelevant  because  we  are  only  interested  in  using  (6.2) 
for  large  values  of  s.  The  upper  root  is  so  =  v£  —  p,  and  by  varying  v  over  a 
suitable  range  we  could  examine  the  behavior  of  Ks,  provided  it  can  be  assumed 
to  vary  slowly  (see  table  I),  after  substituting  the  known  approximation  (5.9) 
in  (6.1)  and  (6.2)  to  cover  the  lower  range  of  s. 

We  shall  adopt  an  approach  which  is  based  on  this  idea  but  is  rather  more 
direct  for  the  present  problem.  Let 

(6.5)  Ht  —  1  +  s/(p  +  s)  +  (  —  )s(s/p)aLs. 

We  have  seen  that  the  effect  of  La  can  be  ignored  at  least  for  s  <  pe/4.  On 
substituting  (6.5)  in  (4.5),  we  get 


(i)fcr.)'thr 


+  S  Lt 

8=0 


=  A  -\-  B, 

say.  The  second  term  B  can  be  expressed  as  2Z  LSUS{  1).  As  before,  put  z  =  s/£ 
and  U8(  1)  =  U(z,  1  )dz.  The  upper  root  of  (6.4)  is  s0  =  £  —  p  and  provided  this  is 
far  from  zero,  U(z,  1)  will  have  an  isolated  peak  at  z0  =  1  —  p/£  near  which 

(6.7)  U{z,  1)  ~  [%/2irz0(l  -  z0)]l/2  exp{  — ^[1  +  c'(z„)]2(2  -  Zo)2/2z0(l  -  z0)}. 
Assuming  that  Ls  =  L(z )  varies  slowly  near  z0,  we  can  infer  that 

(6.8)  B  ~  L(z0)/[1  +  c'(2o)]  =  pLi-p/(S  -  p). 
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Consider  first  the  ease  where  £  is  much  greater  than  p  which  is  itself  large. 
Under  these  conditions  there  is  a  remarkably  simple  approximation  for  p(.r|£, 
r /,  p)  in  the  region  of  large  epidemics.  The  terms  in  A  die  away  rapidly  and 
their  sum  approximates  to 

(G.9)  A  ~  £  (  —  £e~t/p)Vs!  =  exp  (  —  £e~*/p). 

s  =0 

Hence, 

(6.10)  pL*_p/(£  -  p)  ~  -exp  (  — £e_*/p). 

The  formulae  for  p(0|£,  77,  p)  corresponding  to  (6.6)  is 

(cm  P(oif,  „  p)  ~|o  (1  +  jfj)  (-)•  (*) 


s  =  0 


t—fi+V 


Suppose  that  77  is  small.  Compared  with  (6.6),  to  the  order  of  approximation 
considered  the  effect  of  the  extra  factor  [p/(p  +  s)]v  is  to  leave  the  first  term 
unaltered  and  to  multiply  the  second  term  by  [p/(p  +  s,,)]7’  =  (p/£)’1.  It  follows 
that 

(6.12)  p(0|£,  77,  p)  ~  [1  ~  (p/f)"]  exp  (-£<r*/p). 

Then  from  (4.3)  we  obtain,  for  small  values  of  x, 

(6.13)  p{x |£,  77,  p)  ~  [1  —  (p/£)’]  ^  exp  (-£e~*/p). 


In  other  words,  when  the  threshold  is  large  but  the  population  size  is  much 
larger,  the  distribution  of  the  number  remaining  uninfected  in  a  large  epidemic  has 
approximately  the  Poisson  form  uith  the  deterministic  mean  £e~  */p. 

This  could  perhaps  have  been  conjectured  from  (4.3)  on  the  plausible  sup¬ 
position  that  p(0|£  —  x,  77,  p  +  x)  changes  slowly  with  x  when  x  is  small  and  £ 
is  far  above  the  threshold  p.  It  is  a  good  approximation  for  large  values  of  £/p, 
but  otherwise  it  is  rather  a  crude  fit.  One  feels  that  there  must  be  a  direct 
argument  in  terms  of  the  epidemic  process  itself  to  explain  this  Poissonlike 
behavior,  just  as  the  approximating  birth  and  death  process  explains  the  be¬ 
havior  for  small  epidemics. 

At  the  other  end  of  the  distribution  near  x  =  £,  (5.14)  or  (5.15)  provides  a 
good  approximation.  Notice  that  (5.15)  can  also  be  expressed  as 


(6.14) 


(p/I)’ 


rjp 


(p  +  £)2u,+” 


(2w  +  77  —  1)! 
w\(w  +  77) ! 


in  accordance  with  the  point  of  view  adopted  by  Kendall  and  W  hittle,  (p/£)’’ 
being  the  probability  that  a  large  epidemic  will  not  develop. 
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7.  A  more  refined  approximation 

When  £  is  not  much  larger  than  p  the  argument  of  the  previous  section  breaks 
down  principally  because  the  peak  of  U(z,  1)  at  1  —  p/£  is  no  longer  sufficiently 
isolated  from  the  one  at  z  =  0,  and  the  overlap  begins  to  matter  near  the  tail  of 
the  distribution  when  £  —  x  approaches  p  +  x.  Provided  this  does  not  happen, 
we  can  still  improve  the  previous  approximation  considerably  by  looking  again 
at  the  exact  expression  for  the  first  term  of  (6.6)  which  is,  from  (5.10), 

(7.1)  A  =  . ~yj  f  u^~1(2  -  m/£)(  1  -  e~ul pye-u  du 

—  t —  -  -rr-.  /  vf-1(2  —  t>)(l  —  e~vil  p)(e~(v  dv. 

(£  ~  1)  •  Jo 

TABLE  II 

Distribution  of  the  Final  Number  x  of  Uninfected 
p  is  exact;  p\  is  from  (7.4);  pi  is  from  (6.13); 
only  the  lower  end  of  the  distribution  is  tabulated. 


o 

o 

o 

i-H 

II 

’7=1 

p  =  200 

£  =  100 

rj  =  1 

P  =  25 

Exact 

Best  Fit 

Poisson 

Exact 

Best  Fit 

Poisson 

X 

P 

Pi 

P2 

P 

Pi 

P* 

msm 

BH 

— 

.13363 

.13341 

.12012 

1 

.19039 

.18991 

.22001 

2 

#9 

■SH 

Sh@b 

.16379 

.16321 

.20148 

3 

.04889 

.04882 

.04834 

.11187 

.11134 

.12303 

4 

.07646 

.07636 

.08142 

.06741 

.06700 

.05633 

5 

.09857 

.09843 

.10972 

.03781 

.03752 

.02063 

6 

.10906 

.10890 

.12322 

.02036 

.02017 

.00630 

7 

.10649 

.10632 

.11860 

.01073 

.01061 

.00165 

8 

.09364 

.09349 

.09989 

.00561 

.00553 

.00031 

9 

.07532 

.07520 

.07479 

.00293 

.00288 

.00008 

10 

.05609 

.05599 

.05039 

.00154 

.00151 

.00001 

11 

.03905 

.03898 

.03087 

.00082 

.00080 

— 

12 

.02562 

.02557 

.01733 

.00044 

.00043 

13 

.01595 

.01592 

.00898 

.00024 

.00023 

14 

.00947 

.00945 

.00432 

.00013 

.00013 

15 

.00539 

.00538 

.00194 

.00008 

.00007 

16 

.00296 

.00295 

.00082 

.00003 

.00004 

17 

.00157 

.00156 

.00032 

.00002 

.00002 

18 

.00081 

.00080 

.00012 

.00001 

.00001 

19 

.00040 

.00040 

.00004 

.00001 

.00001 

20 

.00020 

.00020 

.00001 

— 

.00001 

21 

.00009 

.00009 

.00001 

— 

22 

.00004 

.00004 

— 

23 

.00002 

.00002 

— 

24 

.00001 

.00001 

— 
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Applying  the  Laplace  method  of  approximation  to  the  integral,  we  get 


(7.2) 

where 


(2  —  v0)vie^1  ro)(l  —  e  Mf/p)* 

{1  +  (t*  -  1)[(1  +  £/p> o  ~  1]}1/2 


=  f(v o,  £,  P ), 


(7.3)  e*/'-  1  =  (vot/p)/(v o  -  1). 

(When  £/p  is  large,  v0  ~  1  and  this  reduces  to  (6.9)).  Since  y  is  small,  the  effect 
on  A  of  the  extra  factor  [p/(p  +  s)]”  in  (6.11)  is  to  multiply  A  by  vl  to  the  same 
order  of  approximation.  By  the  previous  argument,  we  then  get  for  large  epi¬ 
demics 


(7.4)  p(*|{,  4,  p)  ~  (^)S '+”  (*)  [p3  -  (prl)"]/**.  £  -  *.  P  +  *). 


where  v0  satisfies  (7.3)  with  £  —  x  for  £  and  p  +  x  for  p.  Table  II  shows  (7.4) 
to  be  a  remarkably  good  fit  even  for  values  as  low  as  p  =  25,  £  =  100.  When 
p  =  50,  £  =  100  the  fit  is  still  found  to  be  good  for  small  values  of  x  but  it 
begins  to  deteriorate  in  the  tail  of  the  distribution  because  of  the  overlap  effect 
mentioned.  There  seems  to  be  no  simple  way  of  allowing  for  this  with  the 
present  technique. 

From  the  practical  point  of  view,  agreement  to  such  a  high  order  of  accuracy 
is  not  particularly  important  because  the  underlying  model  is  itself  very  much 
idealized,  and  exact  computations  can  in  any  case  be  carried  out  by  computer 
on  the  original  backward  or  forward  equations.  But  it  does  give  one  considerable 
confidence  in  the  method  used  to  arrive  at  the  approximations. 


I  am  much  indebted  to  Mr.  II.  L.  Holder  and  Dr.  V.  D.  Barnett  for  calcula¬ 
tions  on  our  English  Electric  KDF.9.  computer  which  have  sustained  me 
throughout  the  work,  and  in  particular  for  the  computations  of  tables  I  and  II. 
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